XeF2 1 - 300eV Multi-orb processing

PH

  • v2b 27/10/20

    • Updating plot styles to HV.
    • Added branching ratios & more literature data.
    • Added ion yield data.
  • v1 27/03/20

    • Running manually on AntonJr, ePSprc-v1.2 env.
    • Code adapted from recent multi-E chuncked job notebooks.

For previous XeF2 processing, see here.

Setup

In [1]:
!hostname
Stimpy
In [2]:
!conda env list
# conda environments:
#
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\dataVis3D
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\dataVis3D-yt
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSdev
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ePSpkgTest2
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\fibre-sim
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\ipykernel_py2
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\pkgTest
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\pypi-test
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\seabreeze-bk
                         C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\envs\test
base                     C:\ProgramData\Anaconda3
ePSdev                   C:\Users\femtolab\.conda\envs\ePSdev
epsInstallTest           C:\Users\femtolab\.conda\envs\epsInstallTest
epsdev                *  C:\Users\femtolab\.conda\envs\epsdev
seabreeze                C:\Users\femtolab\AppData\Local\conda\conda\envs\seabreeze
tensorflow               C:\Users\femtolab\AppData\Local\conda\conda\envs\tensorflow

In [3]:
import sys
import os
from pathlib import Path
import numpy as np
# import epsproc as ep
import xarray as xr

import matplotlib.pyplot as plt

from datetime import datetime as dt
timeString = dt.now()
In [4]:
# For module testing, include path to module here, otherwise use global installation
if sys.platform == "win32":
    modPath = r'D:\code\github\ePSproc'  # Win test machine
    winFlag = True
else:
    modPath = r'/home/femtolab/github/ePSproc/'  # Linux test machine
    winFlag = False
    
sys.path.append(modPath)
import epsproc as ep

# Plotters
from epsproc.plot import hvPlotters

# Multijob class dev code
from epsproc.classes.multiJob import ePSmultiJob
* pyevtk not found, VTK export not available. 
In [5]:
hvPlotters.setPlotters(width = 700)
# import bokeh
# import holoviews as hv
# hv.extension('bokeh')

Load data

In [6]:
# # Scan for subdirs, based on existing routine in getFiles()

fileBase = r'D:\projects\ePolyScat\xef2\XeF2_highRes_wf'  # Data dir on Stimpy
In [7]:
data = ePSmultiJob(fileBase, verbose = 0)
In [8]:
keys = [4,5,6]  # Set for 4d data only
data.scanFiles(keys = keys)
data.jobsSummary()
Found 3 directories, with 80 files.

*** Job orb21 details
Key: orb21
Dir D:\projects\ePolyScat\xef2\XeF2_highRes_wf\orb21_A1G, 30 files.
{   'batch': 'ePS XeF2_2020, batch XeF2_highRes_wf, orbital orb21_A1G',
    'event': 'orb 21 ionization (Xe 4d, A1G/SG), sph grid. Inputs based on '
             'original 2019 calcs, now with chunking for higher E resolution.',
    'orbE': -76.581003676086,
    'orbLabel': 'Xe 4d, A1G/SG'}

*** Job orb22 details
Key: orb22
Dir D:\projects\ePolyScat\xef2\XeF2_highRes_wf\orb22_E1G, 30 files.
{   'batch': 'ePS XeF2_2020, batch XeF2_highRes_wf, orbital orb22_E1G',
    'event': 'orb 22/23 ionization (Xe 4d, E1G/PG), sph grid. Inputs based on '
             'original 2019 calcs, now with chunking for higher E resolution.',
    'orbE': -76.27895729126399,
    'orbLabel': 'Xe 4d, E1G/PG'}

*** Job orb24 details
Key: orb24
Dir D:\projects\ePolyScat\xef2\XeF2_highRes_wf\orb24_E2G, 20 files.
{   'batch': 'ePS XeF2_2020, batch XeF2_highRes_wf, orbital orb24_E2G',
    'event': 'orb 24/25 ionization (Xe 4d, E2G/DG), sph grid. Inputs based on '
             'original 2019 calcs, now with chunking for higher E resolution.',
    'orbE': -75.67486452162,
    'orbLabel': 'Xe 4d, E2G/DG'}

System properties

In [9]:
data.molSummary()
*** Molecular structure
*** Molecular orbital list (from ePS output file)
EH = Energy (Hartrees), E = Energy (eV), NOrbGrp, OrbGrp, GrpDegen = degeneracies and corresponding orbital numbering by group in ePS, NormInt = single centre expansion convergence (should be ~1.0).
props Sym EH Occ E NOrbGrp OrbGrp GrpDegen NormInt
orb
1 SG -1276.2548 2.0 -34728.662023 1.0 1.0 1.0 1.000000
2 SG -202.4761 2.0 -5509.655317 1.0 2.0 1.0 1.000001
3 SU -181.5792 2.0 -4941.021704 1.0 3.0 1.0 1.000000
4 PU -181.5770 2.0 -4940.961839 1.0 4.0 2.0 1.000000
5 PU -181.5770 2.0 -4940.961839 2.0 4.0 2.0 1.000000
6 SG -43.1198 2.0 -1173.349523 1.0 5.0 1.0 1.000001
7 SU -36.1973 2.0 -984.978703 1.0 6.0 1.0 1.000000
8 PU -36.1881 2.0 -984.728358 1.0 7.0 2.0 1.000001
9 PU -36.1881 2.0 -984.728358 2.0 7.0 2.0 1.000001
10 SG -26.3038 2.0 -715.762856 1.0 8.0 1.0 0.819982
11 SU -26.3038 2.0 -715.762856 1.0 9.0 1.0 0.805992
12 SG -25.8711 2.0 -703.988489 1.0 10.0 1.0 1.000000
13 PG -25.8673 2.0 -703.885086 1.0 11.0 2.0 1.000001
14 PG -25.8673 2.0 -703.885086 2.0 11.0 2.0 1.000001
15 DG -25.8578 2.0 -703.626577 1.0 12.0 2.0 1.000000
16 DG -25.8578 2.0 -703.626577 2.0 12.0 2.0 1.000000
17 SG -8.5562 2.0 -232.826061 1.0 13.0 1.0 0.999998
18 SU -6.2727 2.0 -170.688861 1.0 14.0 1.0 1.000000
19 PU -6.2547 2.0 -170.199056 1.0 15.0 2.0 0.999999
20 PU -6.2547 2.0 -170.199056 2.0 15.0 2.0 0.999999
21 SG -2.8143 2.0 -76.581004 1.0 16.0 1.0 0.999993
22 PG -2.8032 2.0 -76.278957 1.0 17.0 2.0 1.000000
23 PG -2.8032 2.0 -76.278957 2.0 17.0 2.0 1.000000
24 DG -2.7810 2.0 -75.674865 1.0 18.0 2.0 1.000000
25 DG -2.7810 2.0 -75.674865 2.0 18.0 2.0 1.000000
26 SG -1.5681 2.0 -42.670174 1.0 19.0 1.0 0.986091
27 SU -1.5619 2.0 -42.501464 1.0 20.0 1.0 0.984049
28 SG -1.0877 2.0 -29.597825 1.0 21.0 1.0 0.998577
29 SU -0.7372 2.0 -20.060234 1.0 22.0 1.0 0.999021
30 PU -0.6738 2.0 -18.335032 1.0 23.0 2.0 0.998286
31 PU -0.6738 2.0 -18.335032 2.0 23.0 2.0 0.998286
32 PG -0.6426 2.0 -17.486037 1.0 24.0 2.0 0.998036
33 PG -0.6426 2.0 -17.486037 2.0 24.0 2.0 0.998036
34 SG -0.5518 2.0 -15.015243 1.0 25.0 1.0 0.999527
35 PU -0.4994 2.0 -13.589366 1.0 26.0 2.0 0.999297
36 PU -0.4994 2.0 -13.589366 2.0 26.0 2.0 0.999297
*** Warning: some orbital convergences outside single-center expansion convergence tolerance (0.01):
[[10.          0.81998167]
 [11.          0.80599244]
 [26.          0.9860914 ]
 [27.          0.98404896]]

Orbitals 21 - 25 are the Xe(4d) core levels, with $A_{1g}$(SG), $E_{1g}$(PG) and $E_{2g}$(DG) components, where the classifcations are in full $D_{\infty h}$ and ePolyScat labels in parenthesis. The electronic structure was computed in Gamess (HF/SPK-QZP-rel), and converged to a bond-length of 1.9373A, and 4d energies as shown in the table above (~75 - 76eV). (See Molecular orbitals below for plots of the orbitals.)

Below: literature values. Comparing orb 21 to the lower 4d SO binding energy of 70.34 eV, the current computational results are off from the experimental binding energies by ~6.2 eV.

Source: Bancroft, G. M., P. A. Malmquist, S. Svensson, E. Basilier, U. Gelius, and K. Siegbahn. “Gas-Phase ESCA Studies of Valence and Core Levels in Xenon Difluoride and Xenon Tetrafluoride.” Inorganic Chemistry 17, no. 6 (June 1978): 1595–99. https://doi.org/10.1021/ic50184a040.

XeF2 lit values

These values also compare well with the current experimental results - shown below from the file XeF2_4d_LV_115eV_binding_energy_extractions_with_atomic_Xe.png

XeF2 binding energies

Plot cross-sections and betas

These are from ePolyScat's getCro function, and are LF (unaligned ensemble) results. This provides a good, if general, overview. More detailed results - per orbital and symmetry - are given in the addendum. The structure/oscillations in some channels at higher energies remain to be explored/explained!

Overview (all symmetries, length gauge)

In [10]:
# Comparitive plot over datasets (all symmetries only)
Etype = 'Ehv'  # Set for Eke or Ehv energy scale
Erange=[50, 300]  # Plot range (full range if not passed to function below)
data.plotGetCroComp(pType='SIGMA', Etype = Etype, Erange = Erange)
In [11]:
# Comparative plot over datasets (all symmetries only)
data.plotGetCroComp(pType='BETA', Etype=Etype, Erange=Erange)

Branching ratios

In [12]:
# Stack XS data to new data structure
# NOTE: this is not quite correct, since it forces all data to same Ehv axis, but OK for a quick hack, and easy to pull out branching ratios
dsXS = xr.Dataset()  # Set blank dataset, this is easier for stacking, probably
lText = []

for key in data.data.keys():
    subset = data.data[key]['XS'].sel({'Sym':'All', 'XC':'SIGMA'})  # Set XS data, all syms only
    
    # NOTE currently missing full dataset resolution for orb24, so try interp. (2eV not 1eV step size)
    # Note dropna to ensure no NaNs, see http://xarray.pydata.org/en/latest/interpolation.html#interpolating-arrays-with-nan
    if key == 'orb24':
        subset = subset.dropna('Eke').interp(Eke = dsXS.Eke, method = 'cubic')
    
#     dsXS[key] = subset.copy()  # Set .copy() for safety here!
    dsXS[data.data[key]['jobNotes']['orbLabel']]  = subset.copy()  # Set .copy() for safety here!
    

# Convert to Xarray & normalise
dsXS = dsXS.to_array().rename({'variable':'channel'}).squeeze()
dsXS['total'] = dsXS.sum('channel')  # Sum over channels

# Normalise...
dsXS = dsXS/dsXS['total']
In [13]:
# Plot
Erange = [75, 250]
dsXS.swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}, Type='L').plot.line(x='Ehv');

These are broadly similar over the different gauges. There's quite a bit of low-E structure, and things smooth out at higher-E, with ~2:1 ratio between SG:PG or DG channels.

TODO: detailed comparison with expt. Below are the results including SO states, which either need to be summed over (expt) or simulated (theory) for direct comparison! It looks like the $\Sigma:\Pi$ or $\Delta$ ratios are somewhat similar as is, maybe, some of the low-E structure.

Experimental branching ratios

Experimental data

  • Set 1: RF processed results, 02/10/20
  • Set 2: DH ion yields, 27/10/20
In [14]:
# Load experimental data
dataPathExpt = Path(r'D:\projects\XeF2_Soleil_2019\RF_data_analysis_021020')

# Set data attribs in dict similar to ePS results structure
dataFiles = {}
dataFiles['SIGMA'] = {'Data':'XS', 'File':r'XeF2_Xe4d_4d32_4d52_cross_sec_all_photon_energies_02102020.dat', 
                      'States':['4d3/2', '4d5/2']}
dataFiles['BETA'] = {'Data':'Beta', 'File':r'XeF2_Xe4d_beta_all_photon_energies_02102020.dat', 
                     'States':['\pi1/2 (4d3/2)', '\delta3/2 (4d3/2)', '\sigma1/2 (4d5/2)', '\pi3/2 (4d5/2)', '\delta5/2 (4d5/2)']}

# Update with ion data
dataFiles2 = {}
dataFiles2['ION-low'] = {'Data':'XS', 'File':r'XeF2_ion_yield_low_energy_cal.txt'}
dataFiles2['ION-high'] = {'Data':'XS', 'File':r'XeF2_ion_yield_high_energy_cal.txt'}
                         
In [15]:
# Read data files and convert to Xarray
# 27/10/20 added quick hack to set 2nd array for ion yield data
dataList = []
dataList2 = []
for key in dataFiles:
    dataIn = np.loadtxt(dataPathExpt/dataFiles[key]['File'])
    
    # Convert to Xarray
    dataXr = xr.DataArray(dataIn[:,1:], dims=('Ehv','State'), coords={'Ehv':dataIn[:,0], 'State':dataFiles[key]['States'][0:dataIn.shape[1]-1]})
    dataXr.attrs = dataFiles[key]
    dataXr.attrs['dataPath'] = dataPathExpt
    dataXr.name = key
    dataList.append(dataXr)

# Stack to Xarray
dataExpt = xr.concat(dataList, "XC")
dataExpt['XC'] = list(dataFiles.keys())

for key in dataFiles2:
    dataIn = np.loadtxt(dataPathExpt/dataFiles2[key]['File'])
    
    # Convert to Xarray
    dataXr = xr.DataArray(dataIn[:,1].squeeze(), dims=('Ehv'), coords={'Ehv':dataIn[:,0]})  # Only 1D datasets in this case
    dataXr.attrs = dataFiles2[key]
    dataXr.attrs['dataPath'] = dataPathExpt
    dataXr.name = key
    dataList2.append(dataXr)

dataExpt2 = xr.concat(dataList2, "XC")
dataExpt2['XC'] = list(dataFiles2.keys())
In [16]:
# Quick plot
marker = 'x'

for item in dataExpt:
    plt.figure() # Force new figure
    item.dropna('State').plot.line(x='Ehv', marker=marker, ls=':')  # Include dropna here to remove empty state dims

    # Fix axis labels
    if item.XC == 'SIGMA':
        plt.ylabel('XS/Mb')
    else:
        plt.ylabel(r"$\beta_{LM}$")

Note dataset labels here - this dataset has XS summed over spin-orbit states (labelled 4d3/2 and 4d5/2), and betas for spin-orbit states.

In [17]:
# Quick plot - Ion yields
marker = 'x'

for item in dataExpt2:
    plt.figure() # Force new figure
    item.plot.line(x='Ehv')  # Include dropna here to remove empty state dims

    # Fix axis labels
    plt.ylabel("Yield")

Comparison plots...

In [18]:
# Compare with computational results

pType = 'SIGMA'
Erange = [50, 250]
pltObj, lText = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, returnHandles = True)

# Add expt data - cross-secions
scale = dataExpt.sel({'XC':pType}).max() / data.dataSets['orb21']['XS'].sel({'XC':pType, 'Type':'L'}).max()  # Set scaling to match computational data
(dataExpt.sel({'XC':pType})/scale).dropna('State').plot.line(x='Ehv', marker=marker, ls=':');

# Add expt data - ion yields
ionData = 'ION-high'
scale = dataExpt2.sel({'XC':ionData}).max() / data.dataSets['orb21']['XS'].sel({'XC':pType, 'Type':'L'}).max()  # Set scaling to match computational data
(dataExpt2.sel({'XC':ionData})/scale).plot.line(x='Ehv');

# Manual legend fix
lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
lText.extend([ionData])
plt.legend(lText)

# Fix axis labels
if pType == 'SIGMA':
    plt.ylabel('XS/Mb')
    plt.title('XS (normalised to computational results)')
else:
    plt.ylabel(r"$\beta_{LM}$")
    plt.title(r"$\beta_{LM}$")

EDIT 27/10/20: now with ion yield data, which should be taken as the reliable experimental XS result. This doesn't show structure, but is broadly similar to the photoelectron measurements. Structure in photoelectron results to be ignored! Some of the previous comments still apply.


Original notes 13/10/20:

This is interesting... the structures in the computational XS look very similar to the experimental results, but on a much smaller energy scale. I suspect there could be many reasons for this in general (both computational and experimental)... but given that the betas match well (see below), this possibly suggests that this is a signature of multi-electron effects in the bound state.

Some reasons this might be the case:

  • The ePS computations are single active electron type calculations, hence do not include any multi-electron/state-mixing effects.
  • The good agreement with the betas suggests that the scattering potential, hence continuum wavefunction and angular overlap integrals, are, however, well modelled in the computations.
  • This, to me, points to the radial overlap integrals in the computations being quite different to reality, hence anything which strongly influences the radial wavefunction of the bound and/or continuum state. (As a side-note, the bond length used in the calculations could also have an effect here, this was 1.937A, pretty close to the lit value of 1.988A; any other state-mixing, e.g. vibronic couplings, could also play a role of course, although these might generally be expected to affect the angular momentum composition too.)
  • Given that the 4d states are mixed, and we also expect multi-electron effects, the bound state seems the most likely candidate for this.

Any suggestions/comments here would be appreciated, since this is rather speculative, and I don't have a good overview of the spectroscopy and general expectations here at all.

In [19]:
# Beta comparison plot over orbs
pType = 'BETA'
pltObj, lText = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, returnHandles = True)

# Add expt data
dataExpt.sel({'XC':pType}).dropna('State').plot.line(x='Ehv', marker=marker, ls=':');

# Manual legend fix
lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
plt.legend(lText)

# Fix axis labels
if pType == 'SIGMA':
    plt.ylabel('XS/Mb')
    plt.title('XS (normalised to computational results)')
else:
    plt.ylabel(r"$\beta_{LM}$")
    plt.title(r"$\beta_{LM}$")

As previously, applying an energy shift to the data gives a very nice agreement for the betas...

In [30]:
Eshift = 9.0  # NOTE - this breaks energy selection code later, for some reason - float type and/or tolerance in slice?

# Beta comparison plot over orbs
pType = 'BETA'
pltObj, lText = data.plotGetCroComp(pType=pType, Etype = Etype, Erange = Erange, Eshift = Eshift, returnHandles = True)

# Add expt data
dataExpt.sel({'XC':pType}).dropna('State').plot.line(x='Ehv', marker=marker, ls=':');

# Manual legend fix
lText.extend(dataExpt.sel({'XC':pType}).dropna('State').coords['State'].data)
plt.legend(lText)

# Fix axis labels
if pType == 'SIGMA':
    plt.ylabel('XS/Mb')
else:
    plt.ylabel(r"$\beta_{LM}$")

Addendum

In [21]:
from epsproc.vol import orbPlot

# chemPath = r'/home/femtolab/python/chem/tools/chemlab/chemlab'  # Linux dev machine
chemPath = r'D:\temp\chemlab\chemlab'  # Win dev machine
out = orbPlot.importChemLabQC(chemPath = chemPath)

# Test class
# filename = r"/home/femtolab/ePS/XeF2/electronic_structure/xef2_SPKrATZP_rel_geom.log"  # XeF2 Gamess file  OK
filename = r"D:\projects\ePolyScat\xef2\electronic_structure\xef2_SPKrATZP_rel_geom.log"

mo = orbPlot.molOrbPlotter(chemPath = chemPath, fileIn = filename)
Import OK: Chemlab module <module 'qc.wavefunction' from 'D:\\temp\\chemlab\\chemlab\\qc\\wavefunction.py'>
Import OK: Chemlab module <module 'qc.wavefunction' from 'D:\\temp\\chemlab\\chemlab\\qc\\wavefunction.py'>
Found electronic structure file: D:\projects\ePolyScat\xef2\electronic_structure\xef2_SPKrATZP_rel_geom.log
Read 3 atoms and 160 MOs
*** Grids set OK
StructuredGrid (0x210288586a8)
  N Cells:	117649
  N Points:	125000
  X Bounds:	-2.937e+00, 2.937e+00
  Y Bounds:	-2.937e+00, 2.937e+00
  Z Bounds:	-2.937e+00, 2.937e+00
  Dimensions:	50, 50, 50
  N Arrays:	0

In [22]:
# Set orbs to plot - currently handles list of strs only!
orbs = [21,22,23,24,25]
orbList = [str(item) for item in orbs]
[mo.calcOrb(orbN = orbN) for orbN in orbs];
In [23]:
# Plot orbs as independent figures
[mo.plotOrb(orbN = [item], interactive = False, showAtoms = False) for item in orbList];

This shows orbitals 21 - 25, the full set of 4d components. Note that (22,23) and (24,25) are degenerate pairs, and are treated as one orbital (with occ = 4) in ePolyScat, as noted previously.

Gauge results comparison

This shows a big change in absolute X-section between lenght (L) and velocity (V) gauge calculations (~order of magnitue), but the feature are appoximately the same in both cases (apologies for the poor plots here!).

In [24]:
data.plotGetCro(pType='SIGMA', Etype=Etype, Erange=Erange)
In [25]:
# Branching ratios
dsXS.swap_dims({'Eke':'Ehv'}).sel(**{Etype:slice(Erange[0], Erange[1])}).plot.line(x='Ehv', col='Type');

Beta parameters are pretty much idendical for both gauges, with a broad, shape-resonance type feature at lower photon energies, and lots of small-scale resonance structure to higher photon energies.

In [26]:
data.plotGetCro(pType='BETA', Etype=Etype, Erange=Erange)

Matrix elements

In [28]:
data.lmPlot(dataType = 'matE', Erange = Erange, Etype = Etype, thres = 0.1, logFlag = False, selDims = {'Type':'L'}, fillna = True)
Plotting data XeF2_highRes_wf.orb21_A1G_E1.0_10.0_141.0eV.inp.out, pType=a, thres=0.1, with Seaborn
Plotting data XeF2_highRes_wf.orb22_E1G_E1.0_10.0_141.0eV.inp.out, pType=a, thres=0.1, with Seaborn
Plotting data XeF2_highRes_wf.orb24_E2G_E1.0_10.0_141.0eV.inp.out, pType=a, thres=0.1, with Seaborn

Versions

In [29]:
import scooby
scooby.Report(additional=['epsproc', 'xarray', 'jupyter'])
Out[29]:
Wed Oct 28 09:29:54 2020 Eastern Daylight Time
OS Windows CPU(s) 32 Machine AMD64
Architecture 64bit RAM 63.9 GB Environment Jupyter
Python 3.7.3 (default, Apr 24 2019, 15:29:51) [MSC v.1915 64 bit (AMD64)]
epsproc 1.3.0-dev xarray 0.15.0 jupyter Version unknown
numpy 1.19.2 scipy 1.3.0 IPython 7.12.0
matplotlib 3.3.1 scooby 0.5.6
Intel(R) Math Kernel Library Version 2020.0.0 Product Build 20191125 for Intel(R) 64 architecture applications
In [ ]: